% clear;
% load('Result_Baseline.mat');
%% Distribution of Wealth-to-income ratio
Stat    =   SS_WealthIncome(PP,SS);

%% QW
QW              =   SS.TrProb.UnitTrMat.ReEval*SS.Dist_fzo.QW_fzo;
QW_V            =   (1-PP.XI)*QW;
QW_Vhat         =   PP.XI*QW;

Dist            =   full([QW_V;QW_Vhat]);
%% QQ
ID              =   SS.DistApp.fzo.State(:,2);
ID_z            =   SS.FunApp.ExoState.State(ID,1);
ID_FI           =   SS.FunApp.ExoState.State(ID,2);
ID_RI           =   SS.FunApp.ExoState.State(ID,3);

z               =   PP.ExoState.IdioInc.Node(ID_z);
FI              =   PP.ExoState.Idio_FI.Node(ID_FI);
RI              =   PP.ExoState.Idio_RI.Node(ID_RI);

w               =   RI*SS.w_H + (1-RI)*SS.w_N;
% Income
Income          =   [SS.DistPolicy.V.zl.*w;SS.DistPolicy.Vhat.zl.*w];
% Wealth
Wealth          =   [SS.DistPolicy.V.b;SS.DistPolicy.Vhat.b];
% Wealth-Income
Wealth_Income   =   Wealth./Income;
WI_NonNaN_Idx   =   isfinite(Wealth_Income);

%% Bar plot 
Bins            =   linspace(0,10,11)';
N               =   length(Bins)+1;
Frac            =   zeros(N,1);

for ii=1:N
    if ii==1
        Frac(ii)        =   sum(Dist((Wealth_Income <= Bins(ii)) ));
    elseif ii==N
        Frac(ii)        =   sum(Dist((Wealth_Income > Bins(ii-1)) ));
    else 
        Frac(ii)        =   sum(Dist((Wealth_Income > Bins(ii-1)) & (Wealth_Income <= Bins(ii)) ) );
    end
end

bar([Bins; 11], Frac)

%% Output the data 
TempDS          =   [Wealth_Income,Dist];
TempDS          =   TempDS(Dist>0, :);
save('TableGraphs/WealthIncomeRatioDist.mat','TempDS');
